Answering Questions: Gather Data: Diffusion, 2021-07

# Check requisite packages are installed.
packages <- c(
  "plotly", 
  "dplyr",
  "ReacTran",
  "deSolve"
)
for (pkg in packages) {
  library(pkg, character.only = TRUE)
}

package 㤼㸱plotly㤼㸲 was built under R version 4.0.5Loading required package: ggplot2
package 㤼㸱ggplot2㤼㸲 was built under R version 4.0.3
Attaching package: 㤼㸱plotly㤼㸲

The following object is masked from 㤼㸱package:ggplot2㤼㸲:

    last_plot

The following object is masked from 㤼㸱package:stats㤼㸲:

    filter

The following object is masked from 㤼㸱package:graphics㤼㸲:

    layout

package 㤼㸱dplyr㤼㸲 was built under R version 4.0.4
Attaching package: 㤼㸱dplyr㤼㸲

The following objects are masked from 㤼㸱package:stats㤼㸲:

    filter, lag

The following objects are masked from 㤼㸱package:base㤼㸲:

    intersect, setdiff, setequal, union

package 㤼㸱ReacTran㤼㸲 was built under R version 4.0.5Loading required package: rootSolve
package 㤼㸱rootSolve㤼㸲 was built under R version 4.0.3Loading required package: deSolve
package 㤼㸱deSolve㤼㸲 was built under R version 4.0.3Loading required package: shape
package 㤼㸱shape㤼㸲 was built under R version 4.0.5

Load

Pulling code almost directly from LM1996-NumPoolComScaling-Results-2021-05.Rmd.

dirViking <- c(
  file.path(
    getwd(), "LCAB_LawMorton1996-NumericalPoolCommunityScaling"
  ),
  file.path(
    getwd(), "LCAB_LawMorton1996-NumericalPoolCommunityScaling2"
  )
)
dirVikingResults <- file.path(
  dirViking, c("results-2021-04", "save-2021-05-10") # Latter not 100% yet.
)
resultFormat <- paste0(
  "run-", 
  "%d", # Combination Number, or CombnNum.
  "-", 
  "%s", # Run Seed.
  ".RDS"
)

2021-04 Data

# Copied from LawMorton1996-NumericalPoolCommunityScaling-Calculation.R
set.seed(38427042)

basal <- c(3, 10, 30, 100, 300, 1000)
consumer <- c(3, 10, 30, 100, 300, 1000) * 2
events <- (max(basal) + max(consumer)) * 2
runs <- 100

logBodySize <- c(-2, -1, -1, 1) # Morton and Law 1997 version.
parameters <- c(0.01, 10, 0.5, 0.2, 100, 0.1)

# Need to rerun seedsPrep to get the random number generation right for seedsRun
seedsPrep <- runif(2 * length(basal) * length(consumer)) * 1E8
seedsRun <- runif(runs * length(basal) * length(consumer)) * 1E8

paramFrame <- with(list(
  b = rep(basal, times = length(consumer)),
  c = rep(consumer, each = length(basal)),
  s1 = seedsPrep[1:(length(basal) * length(consumer))],
  s2 = seedsPrep[
    (length(basal) * length(consumer) + 1):(
      2 * length(basal) * length(consumer))
  ],
  sR = seedsRun
), {
  temp <- data.frame(
    CombnNum = 0,
    Basals = b,
    Consumers = c,
    SeedPool = s1,
    SeedMat = s2,
    SeedRuns = "",
    SeedRunsNum = 0,
    EndStates = I(rep(list(""), length(b))),
    EndStatesNum = 0,
    EndStateSizes = I(rep(list(""), length(b))),
    EndStateSizesNum = NA,
    EndStateAssembly = I(rep(list(""), length(b))),
    EndStateAbundance = I(rep(list(""), length(b))),
    Dataset = "2021-04",
    DatasetID = 1,
    stringsAsFactors = FALSE
  )
  for (i in 1:nrow(temp)) {
    seeds <- sR[((i - 1) * runs + 1) : (i * runs)]
    temp$SeedRuns[i] <- toString(seeds) # CSV
    temp$SeedRunsNum[i] <- length(seeds)
  }
  temp$CombnNum <- 1:nrow(temp)
  temp
})

# Note: n + 2 end states. Failure to finish, failure to obtain state, and state.
for (i in 1:nrow(paramFrame)) {
  resultsList <- list(
    "No Run" = 0,
    "No State" = 0
  )
  resultsSize <- list(
    "0" = 0
  )
  resultsAssembly <- list(
    "No Run" = data.frame(),
    "No State" = data.frame()
  )
  seeds <- unlist(strsplit(paramFrame$SeedRuns[i], ', '))
  for (seed in seeds) {
    fileName <- file.path(
      dirVikingResults[paramFrame$DatasetID[i]],
      sprintf(resultFormat, paramFrame$CombnNum[i], seed)
    )
    
    if (file.exists(fileName)) {
      temp <- load(fileName)
      temp <- eval(parse(text = temp)) # Get objects.
      
      if (is.data.frame(temp)) {
        community <- toString(
          temp[[ncol(temp)]][[nrow(temp)]]
        )
        size <- toString(length(temp[[ncol(temp)]][[nrow(temp)]]))
        
        if (community == "") {
          resultsList$`No State` <- resultsList$`No State` + 1
          resultsSize$`0` <- resultsSize$`0` + 1
          
        } else if (community %in% names(resultsList)) {
          resultsList[[community]] <- resultsList[[community]] + 1
          resultsSize[[size]] <- resultsSize[[size]] + 1
          
        } else {
          resultsList[[community]] <- 1
          resultsAssembly[[community]] <- temp
          
          if (size %in% resultsSize) {
            resultsSize[[size]] <- resultsSize[[size]] + 1
          } else {
            resultsSize[[size]] <- 1
          }
        }
      } else {
        resultsList$`No State` <- resultsList$`No State` + 1
        resultsSize$`0` <- resultsSize$`0` + 1
      }
    } else {
      resultsList$`No Run` <- resultsList$`No Run` + 1
      resultsSize$`0` <- resultsSize$`0` + 1
    }
  }
  
  paramFrame$EndStates[[i]] <- resultsList
  paramFrame$EndStatesNum[i] <- length(resultsList) - 2 # ! No State, No Run
  paramFrame$EndStateSizes[[i]] <- resultsSize
  paramFrame$EndStateSizesNum[i] <- length(resultsSize) - 1 # ! 0
  paramFrame$EndStateAssembly[[i]] <- resultsAssembly
}

2021-05 Data

source(
  file.path(getwd(), 
            "LawMorton1996-NumericalPoolCommunityScaling-Settings2.R")
)

oldNrow <- nrow(paramFrame)

paramFrame <- rbind(paramFrame, with(list(
  b = rep(basal, times = length(consumer)),
  c = rep(consumer, each = length(basal)),
  s1 = seedsPrep[1:(length(basal) * length(consumer))],
  s2 = seedsPrep[
    (length(basal) * length(consumer) + 1):(
      2 * length(basal) * length(consumer))
  ],
  sR = seedsRun
), {
  temp <- data.frame(
    CombnNum = 0,
    Basals = b,
    Consumers = c,
    SeedPool = s1,
    SeedMat = s2,
    SeedRuns = "",
    SeedRunsNum = 0,
    EndStates = I(rep(list(""), length(b))),
    EndStatesNum = 0,
    EndStateSizes = I(rep(list(""), length(b))),
    EndStateSizesNum = NA,
    EndStateAssembly = I(rep(list(""), length(b))),
    EndStateAbundance = I(rep(list(""), length(b))),
    Dataset = "2021-05",
    DatasetID = 2,
    stringsAsFactors = FALSE
  )
  for (i in 1:nrow(temp)) {
    seeds <- sR[((i - 1) * runs + 1) : (i * runs)]
    temp$SeedRuns[i] <- toString(seeds) # CSV
    temp$SeedRunsNum[i] <- length(seeds)
  }
  temp$CombnNum <- 1:nrow(temp)
  temp
})
)

# Note: n + 2 end states. Failure to finish, failure to obtain state, and state.
# Modified from above, but with the abundance recorded.
for (i in (oldNrow + 1):nrow(paramFrame)) {
  resultsList <- list(
    "No Run" = 0,
    "No State" = 0
  )
  resultsSize <- list(
    "0" = 0
  )
  resultsAssembly <- list(
    "No Run" = data.frame(),
    "No State" = data.frame()
  )
  resultsAbund <- list(
    "No Run" = "",
    "No State" = ""
  )
  seeds <- unlist(strsplit(paramFrame$SeedRuns[i], ', '))
  for (seed in seeds) {
    fileName <- file.path(
      dirVikingResults[paramFrame$DatasetID[i]],
      sprintf(resultFormat, paramFrame$CombnNum[i], seed)
    )
    
    if (file.exists(fileName)) {
      temp <- load(fileName)
      temp <- eval(parse(text = temp)) # Get objects.
      
      if (is.list(temp) && "Result" %in% names(temp)) {
        
        if (is.data.frame(temp$Result))
          community <- temp$Result$Community[[nrow(temp$Result)]]
        else 
          community <- temp$Result
        
        size <- toString(length(community))
        
        if (community[1] != "") 
          abund <- toString(temp$Abund[community + 1])
        else 
          abund <- ""
        
        community <- toString(community)
        
        if (community == "") {
          resultsList$`No State` <- resultsList$`No State` + 1
          resultsSize$`0` <- resultsSize$`0` + 1
          
        } else if (community %in% names(resultsList)) {
          resultsList[[community]] <- resultsList[[community]] + 1
          resultsSize[[size]] <- resultsSize[[size]] + 1
          
        } else {
          resultsList[[community]] <- 1
          resultsAssembly[[community]] <- temp
          resultsAbund[[community]] <- abund
          
          if (size %in% resultsSize) {
            resultsSize[[size]] <- resultsSize[[size]] + 1
          } else {
            resultsSize[[size]] <- 1
          }
        }
      } else {
        resultsList$`No State` <- resultsList$`No State` + 1
        resultsSize$`0` <- resultsSize$`0` + 1
      }
    } else {
      resultsList$`No Run` <- resultsList$`No Run` + 1
      resultsSize$`0` <- resultsSize$`0` + 1
    }
  }
  
  paramFrame$EndStates[[i]] <- resultsList
  paramFrame$EndStatesNum[i] <- length(resultsList) - 2 # ! No State, No Run
  paramFrame$EndStateSizes[[i]] <- resultsSize
  paramFrame$EndStateSizesNum[i] <- length(resultsSize) - 1 # ! 0
  paramFrame$EndStateAssembly[[i]] <- resultsAssembly
  paramFrame$EndStateAbundance[[i]] <- resultsAbund
}

Test Data

testRowNums <- nrow(paramFrame)
testRowsToAdd <- c(2, 6) # Make sure to put in numerical order!

paramFrame <- with(
  list(
    basal2 = c(5, 10, 15),
    consumer2 = c(20, 40, 60),
    logBodySize = c(-2, -1, -1, 0),
    parameters = c(0.01, 10, 0.5, 0.2, 100, 0.1)
  ),
  {
    set.seed(3680180)
    seedsPrep2 <- runif(2 * length(basal2) * length(consumer2)) * 1E8
    with(list(
      b = rep(basal2, times = length(consumer2)),
      c = rep(consumer2, each = length(basal2)),
      s1 = seedsPrep2[1:(length(basal2) * length(consumer2))],
      s2 = seedsPrep2[
        (length(basal2) * length(consumer2) + 1):(
          2 * length(basal2) * length(consumer2))
      ]
    ), {
      rbind(
        paramFrame,
        data.frame(
          CombnNum = testRowsToAdd,
          Basals = b[testRowsToAdd],
          Consumers = c[testRowsToAdd],
          SeedPool = s1[testRowsToAdd],
          SeedMat = s2[testRowsToAdd],
          SeedRuns = "",
          SeedRunsNum = 0,
          EndStates = I(rep(list(""), length(testRowsToAdd))),
          EndStatesNum = 0,
          EndStateSizes = I(rep(list(""), length(testRowsToAdd))),
          EndStateSizesNum = NA,
          EndStateAssembly = I(rep(list(""), length(testRowsToAdd))),
          EndStateAbundance = I(rep(list(""), length(testRowsToAdd))),
          Dataset = "Test",
          DatasetID = max(paramFrame$DatasetID) + 1,
          stringsAsFactors = FALSE
        )
      )
    }
    )
  }
)

testRowNums <- (testRowNums + 1):nrow(paramFrame)
resultsList <- list(
  list(
    "No Run" = 0,
    "No State" = 0,
    "2, 4, 6, 12, 29" = 1,
    "2, 4, 6, 13, 29" = 1
  ),
  list(
    "No Run" = 0,
    "No State" = 0,
    "8, 10, 12, 14, 15, 16, 39, 43" = 1,
    "8, 12, 14, 15, 16, 38, 39" = 1
  )
)
resultsSize <- list(
  list(
    "0" = 0,
    "5" = 2
  ),
  list(
    "0" = 0,
    "8" = 1,
    "7" = 1
  )
)
resultsAbund <- list(
  list(
    "No Run" = "",
    "No State" = "",
    "2, 4, 6, 12, 29" = "742.88553671712, 80.579233072626, 162.128399850253, 20.2082198699389, 18.8589490510429",
    "2, 4, 6, 13, 29" = "668.664143581837, 119.024146851052, 127.680269383867, 30.657960866033, 13.4844194707944"
  ),
  list(
    "No Run" = "",
    "No State" = "",
    "8, 10, 12, 14, 15, 16, 39, 43" = "20.7665807606491, 32.4461165261454, 80.4033387818895, 817.879722033354, 121.136570782828, 18.0390671088957, 12.3834561177271, 19.9674718543196",
    "8, 12, 14, 15, 16, 38, 39" = "82.592048492812, 138.267379166014, 938.158436379166, 51.8610963745021, 5.03556251837491, 14.1019343145825, 25.9231062711228"
  )
)

for (j in seq_along(testRowNums)) {
  i <- testRowNums[j]
  paramFrame$EndStates[[i]] <- resultsList[[j]]
  paramFrame$EndStatesNum[i] <- length(resultsList[[j]]) - 2 # ! No State, No Run
  paramFrame$EndStateSizes[[i]] <- resultsSize[[j]]
  paramFrame$EndStateSizesNum[i] <- length(resultsSize[[j]]) - 1 # ! 0
  paramFrame$EndStateAbundance[[i]] <- resultsAbund[[j]]
}

Plot

# X, Y, Basal and Consumer.
# Z = Sizes of the Endstates.

plotScalingData <- data.frame(
  CombnNum = rep(paramFrame$CombnNum, paramFrame$EndStatesNum),
  Basals = rep(paramFrame$Basals, paramFrame$EndStatesNum),
  Consumers = rep(paramFrame$Consumers, paramFrame$EndStatesNum),
  Dataset = rep(paramFrame$Dataset, paramFrame$EndStatesNum),
  DatasetID = rep(paramFrame$DatasetID, paramFrame$EndStatesNum)
)

# Communities
comms <- unlist(lapply(paramFrame$EndStates, names))
freqs <- unlist(paramFrame$EndStates)
asmbl <- unlist(paramFrame$EndStateAssembly, recursive = FALSE)
asmbl <- asmbl[comms != "No Run" & comms != "No State"]
freqs <- freqs[comms != "No Run" & comms != "No State"]
comms <- comms[comms != "No Run" & comms != "No State"]

asmbl <- lapply(asmbl, function(d) {
  if (is.null(d)) return(NA)
  if ("Result.Outcome" %in% names(d))
    d %>% dplyr::filter(Result.Outcome != "Type 1 (Failure)" & 
                          Result.Outcome != "Present")
  else
    d$Result %>% dplyr::filter(Outcome != "Type 1 (Failure)" & 
                                 Outcome != "Present")
})

plotScalingData$Communities <- comms
plotScalingData$CommunityFreq <- freqs
plotScalingData$CommunitySeq <- asmbl

# Community Size
temp <- unlist(lapply(strsplit(plotScalingData$Communities, ','), length))
plotScalingData$CommunitySize <- temp

# For usage by the reader.

plotScaling <- plotly::plot_ly(
  plotScalingData,
  x = ~Basals,
  y = ~Consumers,
  z = ~CommunitySize,
  color = ~Dataset,
  colors = c("red", "blue", "black")
)

plotScaling <- plotly::add_markers(plotScaling)

plotScaling <- plotly::layout(
  plotScaling,
  scene = list(
    xaxis = list(type = "log"),
    yaxis = list(type = "log"),
    camera = list(
      eye = list(
        x = -1.25, y = -1.25, z = .05
      )
    )
  )
)

plotScaling

Abundances

# > runif(1) * 1E8
# [1] 82598679
set.seed(82598679)

mats <- list()
poolsall <- list() # name pools used in save data; be careful!

for (i in 1:length(dirViking)) {
  temp <- load(file.path(
    dirViking[i], 
    paste0("LawMorton1996-NumericalPoolCommunityScaling-PoolMats", 
           if (i > 1) i else "", 
           ".RDS")
  ))
  mats[[i]] <- eval(parse(text = temp[1]))
  poolsall[[i]] <- eval(parse(text = temp[2]))
}
pools <- poolsall

# Add in the test datasets.
poolsTemp <- list()
matsTemp <- list()
for (r in testRowNums) {
  testRowRow <- paramFrame[r, ]
  poolsTemp[[testRowRow$CombnNum]] <- with(testRowRow,
    RMTRCode2::LawMorton1996_species(
      Basal = Basals,
      Consumer = Consumers,
      Parameters = c(0.01, 10, 0.5, 0.2, 100, 0.1),
      LogBodySize = c(-2, -1, -1, 0),
      seed = SeedPool
    )
  )
  matsTemp[[testRowRow$CombnNum]] <- with(testRowRow,
    RMTRCode2::LawMorton1996_CommunityMat(
      Pool = poolsTemp[[CombnNum]],
      Parameters = c(0.01, 10, 0.5, 0.2, 100, 0.1),
      seed = SeedMat
    )
  )
}

executing %dopar% sequentially: no parallel backend registered

pools[[i + 1]] <- poolsTemp
mats[[i + 1]] <- matsTemp

oldCandidateData <- load(file.path(getwd(), "candidateDataSoFar.Rdata"))
oldCandidateData <- eval(parse(text = oldCandidateData))

candidateData <- plotScalingData %>% dplyr::group_by(
  CombnNum, Dataset
) %>% dplyr::mutate(
  OtherSteadyStates = dplyr::n() - 1
) %>% dplyr::filter(
  OtherSteadyStates > 0
)
candidateData %>% dplyr::select(-CommunitySeq)

# First, check if it is in the paramFrame.
# Second, check if it is in the saved data from the previous.
# Otherwise, ignore it, we'll figure out what it is and why it is missing later.

candidateData$CommunityAbund <- ""

for (r in 1:nrow(candidateData)) {
  # ID 1:4 are used to identify paramFrame, 5 used to identify abundance
  ID <- candidateData[r, 1:6]
  paramFrameRow <- paramFrame %>% dplyr::filter(
    CombnNum == ID$CombnNum,
    Basals == ID$Basals,
    Consumers == ID$Consumers,
    Dataset == ID$Dataset
  )
  
  if (is.list(paramFrameRow$EndStateAbundance[[1]])) {
    entry <- which(ID$Communities == names(paramFrameRow$EndStateAbundance[[1]]))
    if (length(entry)) {
      candidateData$CommunityAbund[r] <- paramFrameRow$EndStateAbundance[[1]][[entry]]
      next()
    }
  }
  
  if (ID$Dataset == "2021-04") {
    
    oldCandDatRow <- oldCandidateData %>% dplyr::filter(
      CombnNum == ID$CombnNum,
      Basals == ID$Basals,
      Consumers == ID$Consumers,
      Communities == ID$Communities
    )
    
    if (nrow(oldCandDatRow) > 0) {
      if (oldCandDatRow$CommunityAbund != "") {
        candidateData$CommunityAbund[r] <- oldCandDatRow$CommunityAbund
      }
    }
  }
}

for (r in 1:nrow(candidateData)) {
  if (!(candidateData$CommunityAbund[r] == "Failure" |
      candidateData$CommunityAbund[r] == "")) next

  # Random guesses, starting from structured.
  temp <- with(
    candidateData[r, ],
    RMTRCode2::FindSteadyStateFromEstimate(
      Pool = pools[[DatasetID]][[CombnNum]],
      InteractionMatrix = mats[[DatasetID]][[CombnNum]],
      Community = Communities,
      Populations = ifelse(
        pools[[DatasetID]][[CombnNum]]$Type[
          RMTRCode2::CsvRowSplit(Communities)
        ] == "Basal",
        1000, 10)
    )
  )

  if (any(temp) < 1E-4) {
    temp <- "EstimateFailure"
  } else {
    temp <- toString(temp)
  }
  candidateData$CommunityAbund[r] <- temp
}

DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
      taken on this call before reaching TOUT     
In above message, I1 = 5000
 
In above message, R1 = 2843.06
 

an excessive amount of work (> maxsteps ) was done, but integration was not successful - increase maxstepsReturning early. Results are accurate, as far as they go

DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
      taken on this call before reaching TOUT     
In above message, I1 = 5000
 
In above message, R1 = 2876.33
 

an excessive amount of work (> maxsteps ) was done, but integration was not successful - increase maxstepsReturning early. Results are accurate, as far as they go

DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
      taken on this call before reaching TOUT     
In above message, I1 = 5000
 
In above message, R1 = 2870.65
 

an excessive amount of work (> maxsteps ) was done, but integration was not successful - increase maxstepsReturning early. Results are accurate, as far as they go

DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
      taken on this call before reaching TOUT     
In above message, I1 = 5000
 
In above message, R1 = 2873.23
 

an excessive amount of work (> maxsteps ) was done, but integration was not successful - increase maxstepsReturning early. Results are accurate, as far as they go

DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
      taken on this call before reaching TOUT     
In above message, I1 = 5000
 
In above message, R1 = 2871.69
 

an excessive amount of work (> maxsteps ) was done, but integration was not successful - increase maxstepsReturning early. Results are accurate, as far as they gocoercing argument of type 'double' to logicalcoercing argument of type 'double' to logicalcoercing argument of type 'double' to logicalcoercing argument of type 'double' to logicalcoercing argument of type 'double' to logicalcoercing argument of type 'double' to logical

DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
      taken on this call before reaching TOUT     
In above message, I1 = 5000
 
In above message, R1 = 6597.38
 

an excessive amount of work (> maxsteps ) was done, but integration was not successful - increase maxstepsReturning early. Results are accurate, as far as they gocoercing argument of type 'double' to logical

DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
      taken on this call before reaching TOUT     
In above message, I1 = 5000
 
In above message, R1 = 3592.76
 

an excessive amount of work (> maxsteps ) was done, but integration was not successful - increase maxstepsReturning early. Results are accurate, as far as they go

DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
      taken on this call before reaching TOUT     
In above message, I1 = 5000
 
In above message, R1 = 3780.53
 

an excessive amount of work (> maxsteps ) was done, but integration was not successful - increase maxstepsReturning early. Results are accurate, as far as they go

DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
      taken on this call before reaching TOUT     
In above message, I1 = 5000
 
In above message, R1 = 3786.89
 

an excessive amount of work (> maxsteps ) was done, but integration was not successful - increase maxstepsReturning early. Results are accurate, as far as they go

DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
      taken on this call before reaching TOUT     
In above message, I1 = 5000
 
In above message, R1 = 3784.07
 

an excessive amount of work (> maxsteps ) was done, but integration was not successful - increase maxstepsReturning early. Results are accurate, as far as they go

DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
      taken on this call before reaching TOUT     
In above message, I1 = 5000
 
In above message, R1 = 3778.43
 

an excessive amount of work (> maxsteps ) was done, but integration was not successful - increase maxstepsReturning early. Results are accurate, as far as they go

DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
      taken on this call before reaching TOUT     
In above message, I1 = 5000
 
In above message, R1 = 3777.75
 

an excessive amount of work (> maxsteps ) was done, but integration was not successful - increase maxstepsReturning early. Results are accurate, as far as they go

DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
      taken on this call before reaching TOUT     
In above message, I1 = 5000
 
In above message, R1 = 3774.98
 

an excessive amount of work (> maxsteps ) was done, but integration was not successful - increase maxstepsReturning early. Results are accurate, as far as they go

DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
      taken on this call before reaching TOUT     
In above message, I1 = 5000
 
In above message, R1 = 3781.84
 

an excessive amount of work (> maxsteps ) was done, but integration was not successful - increase maxstepsReturning early. Results are accurate, as far as they gocoercing argument of type 'double' to logical

candidateData <- candidateData %>% dplyr::filter(CommunityAbund != "",
                                                 CommunityAbund != "Failure",
                                                 CommunityAbund != "EstimateFailure")

candidateData$CommunityProd <- NA
for (r in 1:nrow(candidateData)) {
  candidateData$CommunityProd[r] <- with(
    candidateData[r, ], 
    RMTRCode2::Productivity(
      Pool = pools[[DatasetID]][[CombnNum]], 
      InteractionMatrix = mats[[DatasetID]][[CombnNum]], 
      Community = Communities, 
      Populations = CommunityAbund
    )
  )
}

Simple Pipe Results

Setup

functionCommunity <- function(
    time, 
    Concentration, 
    parameters, 
    GridThickness,# = Grid,
    Pool,# = preprocessed$redDisPool, 
    Mat,# = preprocessed$redComMat, 
    Area = 1, 
    Diffusion = 1, 
    Verbose = FALSE,
    ExtinctionThreshold = 0, 
    Circular = FALSE,
    nonUniformDiffusion = FALSE
) {
    if(Verbose) print(Concentration)
    
    Concentration[Concentration < ExtinctionThreshold] <- 0
    
    # Columns are per species pipes,
    # Rows are the same point on a pipe but for each species.
    pipeMatrix <- matrix(data = Concentration, ncol = nrow(Pool))
    # Note that matrix automatically fills downward!
    
    # Idea with circular: attach duplicates of first, last evaluation points.
    # Then evaluate as normal. Before finishing, remove the duplicate entries.
    # Order of additions matters.
    if (Circular > 0) {
      pipeMatrix <- rbind(pipeMatrix[-Circular:-1 + 1 +nrow(pipeMatrix), ],
                          pipeMatrix,
                          pipeMatrix[1:Circular, ])
      
      if ("dx" %in% names(GridThickness)) {
        GridThickness <- GridThickness$dx
      }
      GridThickness <- c(GridThickness[-Circular:-1 + 1 + length(GridThickness)],
                         GridThickness,
                         GridThickness[1:Circular])
    }
    
    if(Verbose) print(pipeMatrix)
    
    pipeList <- lapply(1:ncol(pipeMatrix), 
                       FUN = function(i, m) {m[, i]}, 
                       m = pipeMatrix)
    
    # Take the species columns to be in the same order as the pool.
    pipeListTran <- lapply(
      seq_along(pipeList), 
      FUN = function(i, C, A, D, G, M) {
        ReacTran::tran.1D(
          C = C[[i]], A = A,
          D = D * M ^ (
            if (nonUniformDiffusion == "positive") 0.62
            else if (nonUniformDiffusion == "negative") -0.62
            else 0
          ),
          dx = G
        )$dC
      }, 
      C = pipeList, A = Area, D = Diffusion, G = GridThickness, M = Pool$Size)
    
    # Diffusion
    pipeMatrixTran <- matrix(
      unlist(pipeListTran),
      nrow = nrow(pipeMatrix),
      ncol = ncol(pipeMatrix)
    )
    
    if(Verbose) print(pipeMatrixTran)
    
    # Point-wise community dynamics
    pipeMatrixComm <- t(apply(
      pipeMatrix, MARGIN = 1, FUN = function(C, t, P, M) {
      unlist(RMTRCode2::GeneralisedLotkaVolterra(t = t, y = C, 
                                      parms = list(
                                        a = M,
                                        r = P$ReproductionRate
                                      )))
    }, t = time, P = Pool, M = Mat))
    
    if(Verbose) print(pipeMatrixComm)
    
    # Order of deletions does not matter.
    if(Circular) pipeMatrixTran <- pipeMatrixTran[
      c(-1:-Circular, -1:-Circular + 1 + nrow(pipeMatrixTran)), ] 
    if(Circular) pipeMatrixComm <- pipeMatrixComm[
      c(-1:-Circular, -1:-Circular + 1 + nrow(pipeMatrixComm)), ] 
    
    list(pipeMatrixTran + pipeMatrixComm)
  }

pipeFUN <- function(i, dat, pool, mat, 
                    GridL = 1, GridCells = 100,
                    Time = seq(0, 10000, by = 2),
                    nonUniformDiffusion = FALSE) {
  temp <- dat[i, ]
  
  preprocessed <- RMTRCode2::IslandPreprocess(
    Pool = pool, InteractionMatrix = mat,
    Communities = c(
      list(temp$Communities[1]),
      list(temp$Communities[2])
    ),
    Populations = c(
      list(temp$CommunityAbund[1]),
      list(temp$CommunityAbund[2])
    ),
    DispersalPool = 0,
    DispersalIsland = matrix(0, nrow = 2, ncol = 2),
    Tolerance = 0, 
    Verbose = FALSE
  )
  
  # Setup the de system.
  Grid <- ReacTran::setup.grid.1D(L = GridL, N = GridCells)
  
  initialCondition <- with(preprocessed, rbind(
    # Community 1
    matrix(abundance_init[1:length(redCom)], ncol = length(redCom)),
    # Space
    matrix(0, nrow = GridCells - 2, ncol = length(redCom)),
    # Community 2
    matrix(abundance_init[length(redCom) + 1:length(redCom)], ncol = length(redCom))
  ))
  
  # https://www.jstor.org/stable/3071786 figure 1
  # suggests size-dispersal should be 10^0.62
  deSolve::ode.1D(
    y = initialCondition, 
    times = Time,
    func = functionCommunity,
    parms = preprocessed$parameters,
    D = 0.001,
    nspec = length(preprocessed$redCom),
    Verbose = FALSE,
    GridThickness = Grid,
    Pool = preprocessed$redPool, 
    Mat = preprocessed$redComMat, 
    nonUniformDiffusion = nonUniformDiffusion
  )
}

# For each group-dataset,
# For each pair,
# Run Pipe Dynamics,
# Save the result with its pairing
candidateData$TotalID <- paste(candidateData$CombnNum, candidateData$DatasetID)

pipeInteractions <- list()
pipeInteractionsFull <- list()
pipeInteractionsP <- list()
pipeInteractionsPFull <- list()
pipeInteractionsN <- list()
pipeInteractionsNFull <- list()
pipeInteractionsTemp <- list()
pipeInteractionsTempFull <- list()

for (grp in unique(candidateData$TotalID)) {
  candidateDataSubset <- candidateData %>% dplyr::filter(TotalID == grp)
  
  if (nrow(candidateDataSubset) == 1) next()
  
  for (nonUniformDiffusion in c(FALSE, "positive", "negative")) {
    pairingResults <- combn(
      nrow(candidateDataSubset), 2, 
      pipeFUN,
      dat = candidateDataSubset, 
      pool = pools[[
        candidateDataSubset$DatasetID[1]
      ]][[candidateDataSubset$CombnNum[1]]],
      mat = mats[[
        candidateDataSubset$DatasetID[1]
      ]][[candidateDataSubset$CombnNum[1]]],
      simplify = FALSE,
      nonUniformDiffusion = nonUniformDiffusion
    )
    
    pairingResultsAbunds <- lapply(
      pairingResults, function(mat, cells = 100) {
        # Final outcomes, without time.
        mat <- mat[nrow(mat), -1]
        retVal <- list()
        species <- length(mat) / cells
        # The results are copies of the pipe by species
        # We want the end points.
        for (i in c(1, cells)) {
          retVal[[toString(i)]] <- mat[i + cells * ((1:species) - 1)]
        }
        retVal
      }
    )
  
    pipeInteractionsTempFull[[
      toString(nonUniformDiffusion)]][[grp]] <- pairingResults
    pipeInteractionsTemp[[
      toString(nonUniformDiffusion)]][[grp]] <- pairingResultsAbunds
  }
}

pipeInteractionsFull <- pipeInteractionsTempFull[["FALSE"]]
pipeInteractions <- pipeInteractionsTemp[["FALSE"]]
pipeInteractionsPFull <- pipeInteractionsTempFull[["positive"]]
pipeInteractionsP <- pipeInteractionsTemp[["positive"]]
pipeInteractionsNFull <- pipeInteractionsTempFull[["negative"]]
pipeInteractionsN <- pipeInteractionsTemp[["negative"]]

Results Size-Negative

# Format of table should be:
# ID, Community 1, Community 2, Outcomes 1-2, Outcomes 1-0-2
# For outcomes, species presence will be used.

communities <- NULL
totalCommunities <- NULL
for (grp in unique(candidateData$TotalID)) {
  candidateDataSubset <- candidateData %>% dplyr::filter(TotalID == grp)
  
  if (nrow(candidateDataSubset) > 1) {
    newCommunities <- combn(
      candidateDataSubset$Communities, 2, 
    )
    communities <- c(communities, newCommunities)
    totalCommunities <- c( # Labelling is wrong here. Need tC from pairs of nC.
      totalCommunities,
      list(apply(newCommunities, 2, function(coms) {
        toString(sort(unique(RMTRCode2::CsvRowSplit(coms))))
      }))
    )
  }
}

minThresh <- lapply(candidateData$CommunityAbund, function(x) min(RMTRCode2::CsvRowSplit(x)) * 0.1) %>% unlist %>% min

pipeInteractionsNWhich <- unlist(lapply(
  seq_along(pipeInteractionsN), function(i, byID, communities) {
    lapply(byID[[i]], function(system, comms) {
      lapply(system, function(pipeEnd, comms) {
        toString(RMTRCode2::CsvRowSplit(comms)[
          which(pipeEnd > minThresh)
        ])
      }, comms = comms)
    }, comms = communities[[i]])
  }, byID = pipeInteractionsN, communities = totalCommunities
))

pipeInteractionNResults <- data.frame(
  DatasetID = rep(names(pipeInteractionsN), 
                  unlist(lapply(pipeInteractionsN, length))),
  Community1 = communities[seq(from = 1, to = length(communities), by = 2)],
  Community2 = communities[seq(from = 2, to = length(communities), by = 2)],
  Outcome_Pipe1 = pipeInteractionsNWhich[
    seq(from = 1, to = length(pipeInteractionsNWhich), by = 2)],
  Outcome_Pipe2 = pipeInteractionsNWhich[
    seq(from = 1, to = length(pipeInteractionsNWhich), by = 2)]
)

pipeInteractionNResults

pipeInteractionNResults %>% dplyr::mutate(
  C1Invaded = Community1 != Outcome_Pipe1,
  C2Invaded = Community2 != Outcome_Pipe2,
  Stalemate = !C1Invaded & !C2Invaded,
  Hybrid = C1Invaded & C2Invaded,
) %>% dplyr::select(-dplyr::starts_with("Outcome"))

Results Size-Neutral

pipeInteractionsWhich <- unlist(lapply(
  seq_along(pipeInteractions), function(i, byID, communities) {
    lapply(byID[[i]], function(system, comms) {
      lapply(system, function(pipeEnd, comms) {
        toString(RMTRCode2::CsvRowSplit(comms)[
          which(pipeEnd > minThresh)
        ])
      }, comms = comms)
    }, comms = communities[[i]])
  }, byID = pipeInteractions, communities = totalCommunities
))

pipeInteractionResults <- data.frame(
  DatasetID = rep(names(pipeInteractions), 
                  unlist(lapply(pipeInteractions, length))),
  Community1 = communities[seq(from = 1, to = length(communities), by = 2)],
  Community2 = communities[seq(from = 2, to = length(communities), by = 2)],
  Outcome_Pipe1 = pipeInteractionsWhich[
    seq(from = 1, to = length(pipeInteractionsWhich), by = 2)],
  Outcome_Pipe2 = pipeInteractionsWhich[
    seq(from = 1, to = length(pipeInteractionsWhich), by = 2)]
)

pipeInteractionResults

pipeInteractionResults %>% dplyr::mutate(
  C1Invaded = Community1 != Outcome_Pipe1,
  C2Invaded = Community2 != Outcome_Pipe2,
  Stalemate = !C1Invaded & !C2Invaded,
  Hybrid = C1Invaded & C2Invaded,
) %>% dplyr::select(-dplyr::starts_with("Outcome"))

Results Size-Positive

pipeInteractionsPWhich <- unlist(lapply(
  seq_along(pipeInteractionsP), function(i, byID, communities) {
    lapply(byID[[i]], function(system, comms) {
      lapply(system, function(pipeEnd, comms) {
        toString(RMTRCode2::CsvRowSplit(comms)[
          which(pipeEnd > minThresh)
        ])
      }, comms = comms)
    }, comms = communities[[i]])
  }, byID = pipeInteractionsP, communities = totalCommunities
))

pipeInteractionPResults <- data.frame(
  DatasetID = rep(names(pipeInteractionsP), 
                  unlist(lapply(pipeInteractionsP, length))),
  Community1 = communities[seq(from = 1, to = length(communities), by = 2)],
  Community2 = communities[seq(from = 2, to = length(communities), by = 2)],
  Outcome_Pipe1 = pipeInteractionsPWhich[
    seq(from = 1, to = length(pipeInteractionsPWhich), by = 2)],
  Outcome_Pipe2 = pipeInteractionsPWhich[
    seq(from = 1, to = length(pipeInteractionsPWhich), by = 2)]
)

pipeInteractionPResults

pipeInteractionPResults %>% dplyr::mutate(
  C1Invaded = Community1 != Outcome_Pipe1,
  C2Invaded = Community2 != Outcome_Pipe2,
  Stalemate = !C1Invaded & !C2Invaded,
  Hybrid = C1Invaded & C2Invaded,
) %>% dplyr::select(-dplyr::starts_with("Outcome"))

Save workspace

save(
  candidateData,
  pipeInteractions,
  pipeInteractionsWhich,
  pipeInteractionsP,
  pipeInteractionsPWhich,
  pipeInteractionsN,
  pipeInteractionsNWhich,
  mats,
  paramFrame,
  plotScalingData,
  pools,
  file = paste0("LM1996-NumPoolCom-QDatDif",
                "-2021-07.RData")
)

save(
  pipeInteractionsNFull,
  pipeInteractionsFull,
  pipeInteractionsPFull,
  file = paste0("LM1996-NumPoolCom-QDatDif",
                "FULL-2021-07.RData")
)